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ABSTRACT 

The  anthrax  attacks  of  September  18,  2001  refocused 
attention  on  methods  of  asymmetric  warfare  including  in¬ 
troduction  of  contaminants  into  water  distribution  systems. 
To  predict  the  outcome  of  such  an  event,  previous  chemical 
models  of  transport  in  water  distribution  systems  assumed 
that  any  contaminant  would  move  through  the  system  via 
pure  hydraulic  transport,  meaning  that  the  contaminant 
would  not  interact  with  the  pipe  wall.  Such  models  predict 
that  flushing  a  system  with  uncontaminated  water  would 
clear  away  contaminants  fairly  rapidly.  Previous  field  ex¬ 
perience  with  actual  contamination  events  (O’Brien,  2003) 
showed  that  flushing  contaminants  was  quite  slow  and  diffi¬ 
cult  in  practice.  For  example,  in  1980,  a  chlordane  contami¬ 
nation  event  in  Pittsburgh,  PA  required  a  flushing  program 
of  8  months  duration  for  2000  affected  customers.  For  mili¬ 
tary  installations,  such  a  long  delay  represents  an  intolerable 
risk  to  mission  capability  (Hock  et  al.,  2005;  Ginsberg  and 
Hock,  2004).  The  research  described  uses  computational 
chemistry  to  calculate,  ab-initio ,  time  constants  associated 
with  sorption,  desorption,  and  hydrolysis  using  only  selec¬ 
tive  verification  by  empirical  methods.  The  computational 
methods  used  include:  molecular  dynamics  using  NAno- 
scale  Molecular  Dynamics  (NAMD)  to  describe  sorption 
and  desorption,  Turbomol  and  Cosmotherm  to  obtain  corre¬ 
lations  to  molecular  morphology  that  are  less  computation¬ 
ally  expensive  than  NAMD,  and  density  functional  theory 
(DFT)  combined  with  the  polarizable  continuum  model 
(PCM)  in  Gaussian03  to  predict  reaction  rates  of  hydrolysis. 

1.  INTRODUCTION 

In  previous  research,  empirical  methods  showed  that 
uptake  of  contaminant  on  pipe  material  was  a  function  of 
the  contaminant’s  hydrophobicity,  molecular  weight,  and 
number  of  contact  points  with  the  pipe  surface.  Fate  and 
transport  are  also  dependent  on  the  rate  of  hydrolysis  of 
contaminant  in  the  water  distribution  environment,  includ¬ 
ing  interaction  with  chlorination.  This  new  detailed  knowl¬ 
edge  of  fate  and  transport  of  contaminants  in  water  distribu¬ 
tion  systems  indicated  that  it  is  feasible  to  use  computational 


methods  to  directly  calculate  the  time  constants  associated 
with  uptake  and  hydrolysis  using  only  selective  verification 
by  empirical  methods.  All  methods  are  validated  with  labo¬ 
ratory  measurements  (Ginsberg  and  Hock,  2004;  Ginsberg 
et  al.,  2005) 

2.  EXPERIMENTAL  RESULTS 

Previous  laboratory  results  showed  that  uptake  of  con¬ 
taminant  on  pipe  surfaces  is  characterized  by  thermody¬ 
namic  partition,  where  a  fixed  fraction  of  the  contaminant 
available  sorbs  to  the  pipe  material.  The  rate  of  interaction  is 
a  first  order  relaxation  to  an  asymptotic  equilibrium  uptake. 
As  such  it  is  characterized  by  the  asymptotic  fraction  ‘/T  (on 
a  scale  of  0  for  no  uptake,  to  1  for  complete  uptake),  and  an 
exponential  time  constant  T’ .  Laboratory  measurements  are 
being  conducted  on  several  combinations  of  (CB)  contami¬ 
nants  and  pipe  material.  The  results  indicate  that  the  interac¬ 
tion  between  pipe  wall  and  contaminant  can  be  quite  strong. 
As  an  example,  Figure  1  shows  the  asymptotic  uptake  of 
prophos  (a  commercially  available  nematacide)  on  concrete- 
lined  ductile-iron  pipe  with  an  asphaltic  sealer  coat  is  100% 
at  an  exponential  relaxation  time  constant  of  77  hours  (i.e., 
P=1  and  T=77  hours).  As  77  hours  is  comparable  to  the  total 
transport  time  of  water  within  a  distribution  system,  the  ex¬ 
pected  effect  under  field  conditions  would  be  quite  strong. 
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Figure  1 .  Experimental  result  showing  aggressive  uptake 
of  prophos  in  concrete-lined  pipe  with  an  asphaltic  sealer 
coat. 
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3.  MOLECULAR  DYNAMICS  (MD) 
STUDY  OF  FLOW-INDUCED 
SORPTION/DESORPTION. 

To  understand  the  microscopic  mechanism  of  contami¬ 
nant  interaction  with  the  pipe  surfaces  we  performed  MD 
simulations  of  water/pipe/contaminant  systems.  MD  is  a 
computational  method  that  allows  the  microscopic  dynamics 
of  a  molecular  system  to  be  modeled  at  atomic  resolution.  In 
this  method,  a  molecular  system  is  approximated  by  an  en¬ 
semble  of  virtual  atoms  interacting  with  each  other  accord¬ 
ing  to  a  complex  potential  function  -  the  molecular  force 
field  (Allen  and  Tildesley,  1987),  which  has  been  developed 
and  calibrated  to  quantitatively  reproduce  the  physical  prop¬ 
erties  of  a  simulated  system.  The  trajectory  of  each  atom  is 
obtained  by  integrating  Newton’s  equation  of  motion.  The 
potential  energy  function  describing  the  interactions  of  all 
atoms  in  the  system  is  derived  from  quantum-mechanics 
calculations  and  is  further  refined  to  reproduce  physical 
properties  of  the  simulated  system. 

For  our  test  simulations  of  flow-induced  sop- 
tion/desorption  we  used  Dimethyl  Methylphosphonate 
(DMMP)  as  a  prototypical  organophosphate  contaminant. 
As  no  force  field  for  DMMP  consistent  with  our  model  of 
water  and  silica  was  available,  we  have  adapted  the  force 
field  parameters  for  DMMP  based  on  its  rigid-bond  model 
(Vishnyakov  and  Neimark  2004),  additional  parameters 
from  (Barvik  et  al.,  2002)  and  the  CHARMM  force  field 
(MacKerell  1998).  The  obtained  parameters  were  tested  by 
repeating  simulations  using  a  method  reported  in  Vishnya¬ 
kov  and  Neimark  (2004). 

As  a  model  of  typical  surfaces  found  in  a  cement  pipe, 
we  used  an  all-atom  model  of  amorphous  silica  that  accu¬ 
rately  reproduces  atomic-scale  roughness  of  the  surface  and 
its  wetting  properties  (Cruz-Chu  et  al.,  2006).  The  interac¬ 
tions  between  water  and  silicon  dioxide  were  calibrated  by 
measuring  a  water  contact  angle  (WCA),  the  angle  between 
the  tangent  of  a  sessile  water  droplet  and  the  solid  surface 
(Werder  et  al.,  2003).  The  silica  surfaces  were  obtained  us¬ 
ing  numerical  procedures  that  mimic  experimental  anneal¬ 
ing,  varying  the  heating  and  quenching  rates.  The  force  field 
parameters  describing  the  interactions  of  water  and  silica 
were  calibrated  to  reproduce  the  experimental  dependence 
of  the  WCA  on  the  number  of  dangling  atoms  at  the  surface, 
the  latter  being  controlled  by  the  annealing  cycle.  Thus,  our 
model  is  a  close  representation  of  the  real  surface,  providing 
a  good  compromise  between  the  heterogeneity  of  amor¬ 
phous  silica  and  the  magnitude  of  the  intermolecular  forces. 
For  our  studies  of  flow-induced  sorption/desorption,  we 
used  silica  surfaces  of  the  following  two  types:  an  atomi¬ 
cally  flat  hydrophilic  surface  containing  a  large  number  of 
dangling  atoms,  and  an  annealed  hydrophobic  surface  that 
had  no  dangling  atoms  but  had  a  considerable  atomic-scale 
roughness  (Figure  3). 


Figure  2.  Molecular  dynamics  simulations  of  flow- 
induced  sorption/desorption  of  DMMP. 


A  typical  simulation  system,  shown  in  Fig.  2,  contained 
an  aqueous  solution  of  DMMP  placed  between  two  rectan¬ 
gular  slabs  of  silica.  MD  simulations  were  performed  using 
NAMD  (Phillips  et  al.,  2005),  periodic  boundary  conditions, 
and  full  electrostatics.  Other  simulation  conditions  and  pro¬ 
tocols  used  were  as  described  in  (Aksimentiev  et  al.,  2004). 
To  induce  a  flow  of  water  through  the  channels,  each  water 
molecule  was  subject  to  a  constant  force  6F’  directed  along 
the  desired  water  direction.  For  a  system  with  ‘N’  atoms, 
this  creates  a  pressure  difference,  AP  =  NF/A,  where  ‘A  ’  is 
the  cross-section  of  the  channel. 

To  test  if  the  character  of  the  silica  surface  can  influ¬ 
ence  the  uptake  of  the  contaminants,  flow-induced  simula¬ 
tions  were  performed  through  channels  having  hydrophobic 
and  hydrophilic  surfaces.  As  the  mixture  of  DMMP  and 
water  was  pushed  to  flow  through  the  channels,  DMMP 
molecules  were  observed  to  interact  with  the  channel’s  sur¬ 
faces;  some  of  those  interactions  led  to  DMMP  binding. 
After  10  ns  of  MD,  both  systems  reached  a  quasi-stationary 
state.  Figure  3  shows  the  results  of  these  simulations.  In  the 
hydrophilic  channel,  about  70%  of  the  contaminant  ad¬ 
sorbed  to  the  surface;  the  DMMP  molecules  were  observed 
to  bind  to  and  unbind  from  the  surface  in  a  dynamic  equilib¬ 
rium  state.  In  the  hydrophobic  channel,  all  contaminant 
molecules  adsorbed  to  the  surface  after  about  10  ns  simula¬ 
tion.  This  dramatic  influence  of  the  surface  properties  on  the 
update  of  organophosphate  contaminants  suggests  that  the 
latter  may  be  controlled  though  surface  treatment. 

These  simulations  demonstrated  the  feasibility  of  using 
the  MD  method  to  quantitatively  characterize  sorp¬ 
tion/desorption  phenomena  in  a  flow  through  a  channel. 
Future  studies  will  be  directed  toward  computing  the  rates 
of  sorption/desorption,  characterizing  the  influence  of  the 
surface  heterogeneity  on  the  contaminant  binding,  compara¬ 
tive  studies  of  binding  kinetics  of  difference  contaminants 
and  on  extending  this  methodology  to  modeling  adsorption 
of  biomacromolecules. 
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Figure  3.  Surface  properties  affect  the  uptake  of  DPPM.  The  number  of  DMMP  molecules  adsorbed  at 
the  hydrophilic  (left)  and  hydrophobic  (right)  surfaces  is  plotted  versus  simulation  time.  The  snapshots 
at  the  bottom  illustrate  the  two  model  surfaces  used. 


4.  COMPUTATIONAL 
CHARACTERIZATION  USING 
TURBOMOL  AND  COSMOTHERM 

The  theory  of  adsorption  has  far-reaching  applications 
within  the  scientific  and  engineering  community.  One 
branch  of  this  theory  which  has  received  much  interest 
within  the  last  decade  has  been  developed  within  framework 
of  the  conductor-like  screening  (COSMO-RS)  model 
(Klamt  and  Schuurmann,  1993;  Klamt,  1995;  Andzelm  et 
al.,  1995;  MaaBen  et  al.,  1995).  This  theory,  which  com¬ 
bines  aspects  of  quantum  theory,  dielectric  continuum  the¬ 
ory  and  statistical  thermodynamics,  was  originally  devel¬ 
oped  for  fluid  phase  property  prediction.  However,  the  au¬ 
thors  were  quick  to  grasp  the  possibilities  of  extension  to 
partitioning  complex  systems,  and  applications  ranging  from 
intestinal  adsorption  to  soil  sorption  to  adsorption  on  acti¬ 
vated  carbon  rapidly  ensued.  This  portion  of  the  project  fol¬ 
lows  within  the  spirit  of  this  devel¬ 
opment,  and  attempts  to  apply  a 
similar  methodology  to  explain 
partitioning  of  water-borne  con¬ 
taminants  between  the  aqueous 
phase  and  the  pipe  wall.  The  basis 
of  this  model  is  in  the  surface  po¬ 
larization  charge  density  (Figure 
4).  This  quantity  is  known  to  be  a 
very  good  local  measure  of  polar¬ 
ity  and  therefore  correlate  strongly 
with  a  chemical’s  response  to  its 


environment  at  the  molecular  level.  More  general  forms  of 
the  method  require  a  minimum  of  fitted  parameters  and  gen¬ 
eration  of  the  profile  for  all  compounds  involved  in  the  case 
of  a  standard  miscibility  or  partitioning  scenario,  such  as  the 
octanol- water  partition  coefficient.  In  the  case  of  adsorption 
to  a  solid  surface,  less  direct  means  may  be  used,  which 
involve  more  mathematical  fitting  to  sufficient  quantities  of 
data.  The  chemical  (and  to  some  extent  structural)  features 
of  the  pipe  walls  of  interest  for  this  project  are  largely  un¬ 
known,  and  for  this  reason  a  mathematical  fit  will  be  at¬ 
tempted. 

5.  COMPUTATIONAL  PREDICTIONS  OF 
HYDROLYSIS  IN  POTABLE  WATER 

Ab-initio  quantum  chemical  techniques  were  used  to 
investigate  the  degradation  mechanism  of  a  representative 
organophosphate,  Phosphonothioic  acid,  in  water.  Calcula¬ 
tions  using  the  Density  Functional  Theory  (DFT)  approach 


Figure  4.  Polarity  maps  of  the  molecules:  1 -octanol,  DMMP,  and  water.  The 
intense  red  and  blue  areas  are  polar,  green  is  non-polar. 
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and  the  B3LYP  exchange  and  correlations  functionals  were 
performed  to  optimize  the  structure  of  an  organophosphate 
in  the  gas  phase  and  in  water  (Figure  5).  The  Polarizable 
Continuum  Model  (PCM)  approach  available  in  the  Gaus- 
sian03  (Frisch  et  al,  2004)  suite  of  programs  was  used  to 
simulate  the  solvent  water. 


Figure  5.  Optimized  structure  of  an  organophosphate. 


Previous  studies  on  the  hydrolysis  of  phosphinate, 
phosphonate  and  phosphate  molecules  showed  that  good 
accu-racy  was  obtained  using  this  method  and  a  6- 
31+G(2d,2p)  basis  set  (Wright  et  al,  2005),  and  this  basis 
was  used  in  the  ensuing  calculations.  Tight  convergence 
criteria  were  used  throughout  the  study.  The  geometry  ob¬ 
tained  from  the  optimization  in  water  was  used  as  a  starting 
point  for  the  hydrolysis  reaction.  Verification  that  a  transi¬ 
tion  state  was  obtained  in  the  degradation  pathway  was  ob¬ 
tained  by  carefully  tracing  the  reaction  path  in  both  direc¬ 
tions  from  the  transition  structure.  The  predicted  mechanism 
shows  that  an  OH-  group  from  a  water  molecule  approaches 
the  phosphonate  group  opposite  to  the  fluoride  leaving 
group  through  an  SN2-type  reaction.  Competition  exists 
between  a  breaking  of  the  P-S  bond  and  a  breaking  of  the  P- 
OEt  bond  in  the  phosphonate  group,  as  has  been  observed 
by  Daniel  et  al.  (2008).  We  are  currently  working  to  deter¬ 
mine  the  complete  energetics  of  the  system,  and  to  deter¬ 
mine  reaction  rates  for  the  competing  cleavage  reactions. 

6.  DISCUSSION 

The  computational  methods  summarized  here  represent 
several  advances  in  computational  chemistry  as  applied  to 
infrastructure:  1.)  These  computational  tools  are  inexpen¬ 
sive  compared  to  using  purely  empirical  methods  to  obtain 
the  same  results.  2.)  Although  the  theory  behind  scaling 
from  nano- scale  interactions  between  contaminants  and  pipe 
surfaces  has  been,  in  principle,  settled  by  the  statistical  view 
of  thermal  physics,  the  ability  to  make  the  connection  be¬ 
tween  forces  measured  in  molecular  computer  simulations 
and  the  Gibbs  free  energy  is  relatively  recent,  (ex.  Roux, 
2002)  3.)  Once  the  results  are  understood  for  the  computa¬ 
tionally  intensive  studies  using  molecular  dynamics  and 
density  functional  theory,  correlation  studies  may  show  how 
to  obtain  the  same  results  with  tools  like  Turbomol  /  COS- 
MOtherm  using  more  modest  computer  resources. 


This  research  also  establishes  a  scientific  basis  for  sub¬ 
sequent  technical  improvements  in  military  and  civilian 
simulation  capabilities.  The  time  constants  and  partition 
coefficients  will  be  embedded  in  current  hydraulic  simula¬ 
tions  tools  (models)  such  as  EPANet  or  Aquis.  It  represents 
the  first  research  to  predict  the  strong  interactions  between 
chemical  contaminants  and  commonly  used  pipe  material. 
The  improved  fidelity  of  the  models  underlying  water  dis¬ 
tribution  simulations  will  provide  the  Army  with  improved 
tools  for  vulnerability  assessment,  emergency  response 
planning,  and  infrastructure  planning.  Immediate  military 
applications  also  include  development  of  risk  assessment 
models  for  force  projection  and  for  critical  fixed  facilities 
(e.g.,  military  hospitals). 

Beyond  its  military  applications,  this  research  effort 
contributes  substantially  to  the  national  welfare.  In  civilian 
applications,  the  improved  simulation  capability  will  assist 
better  formulation  of  capital  improvement  plans,  facility 
protection  methods,  and  emergency  response  strategies  for 
iconic  structures  and  vulnerable  fixed  facilities  (e.g.,  food 
manufacturing,  judicial,  and  governmental  functions).  This 
helps  formulate  immediate  response  strategies  for  water¬ 
borne  attacks  such  as  those  already  attempted  by  terrorist 
groups. 

7.  CONCLUSION 

The  research  being  conducted  shows  that  three  princi¬ 
ple  computational  approaches,  applied  in  a  synergistic  way 
and  with  empirical  validation,  are  being  used  to  determine 
the  fate  and  transport  constants  of  contaminants  in  a  water 
distribution  system.  Sorption  of  contaminants  on  pipe  mate¬ 
rial  can  be  simulated  using  molecular  dynamics.  Hydrolysis 
rates  are  determined  using  a  combination  of  density  func¬ 
tional  theory  and  the  polarizable  continuum  model.  Both 
sets  of  computational  results  may  later  be  extended  by  dem¬ 
onstrating  that  they  correlated  to  less  computationally  inten¬ 
sive  methods  such  as  Turbomol  /  COSMOtherm.  The  fate 
and  transport  constants,  in  turn,  can  provide  new  simulation 
tools  that  can  enhance  vulnerability  assessment  methods, 
emergency  response  planning,  infrastructure  planning,  and 
system  operation. 
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